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JL ' Abstract 

1^ ' Quarkonium production in high-energy proton (deuteron)-nucleus colhsions is inves- 

tigated in the color glass condensate framework. We employ the color evaporation 
model assuming that the quark pair produced from dense small-x gluons in the nu- 
^ I clear target bounds into a quarkonium outside the target. The unintegrated gluon 

distribution at small Bjorken x in the nuclear target is treated with the Balitsky- 
^^ . Kovchegov equation with running coupling corrections. For the gluons in the proton, 

^v^ I we examine two possible descriptions, unintegrated gluon distribution and ordinary 

collinear gluon distribution. We present the transverse momentum spectrum and 
f^ . nuclear modification factor for 3/ip production at RHIC and LHC energies, and 

^^ I those for T(IS') at LHC energy, and discuss the nuclear modification factor and the 

momentum broadening by changing the rapidity and the initial saturation scale. 



H ■ 1 Introduction 



High-energy proton-nucleus (pA) collisions allow us to explore the dense gluon system 
appearing at small values of the Bjorken x in the target nucleus. Such a dense gluon 
system is expected to possess a universal feature of parton saturation, characterized by 
the saturation momentum scale Ql{x), and has been investigated with the Color Glass 
Condensate (CGC) effective theory[l, 2]. In a heavy nucleus with atomic mass number 
A, the saturation phenomenon will become relevant even at moderate x because of its 
larger gluon density by a factor of target thickness A^^^. Very recently p+Pb collisions at 
the center-of-mass energy ^/s = 5.02 TeV have been delivered at the large hadron collider 
(LHC), and new exciting data are being reported such as hadron multiplicity [3] and ridge 
phenomenon[4, 5]. The observed hadron multiplicity and momentum spectrum [3] will 
constrain theoretical models. 

At this high energy, the relevant value of a;i,2 becomes so small that the parton satu- 
ration scale Q^{x2) oc A^''^X2 with A ^ 0.3[6, 7] in the nuclear target will be comparable 
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to or larger than the charm quark mass rric. This suggests coherence effects even in 
charm quark production, and thus we can get useful information of the gluon saturation 
in the nuclear target by studying the energy and rapidity dependences of charm quark 
and quarkonium productions in pA collisions at relativistic heavy ion collider (RHIC) and 
the LHC[8, 9, 10, 11, 12, 13]. 

In heavy ion collision experiments, heavy quark and quarkonium productions[14, 
15] are very valuable probes for quantifying properties of hot and dense matter or a 
quark-gluon plasma transiently created in the events. Quarkonium suppression[16] and 
enhancement [17, 18], and energy loss[19, 20, 21] and collective flow of heavy flavor mesons[22] 
have already been discussed extensively. Here one obviously needs to know the initial 
nuclear effects on their productions for quantitative understanding of hot-medium modi- 
fications. 

In this paper we shall investigate phenomenological implications of the saturation ef- 
fects on the quarkonium production at collider energies by exploiting the color evaporation 
model (CEM)[15]. Systematic study of quarkonium production spectrum from CGC in 
pA collisions will quantitatively improve our understanding of gluon saturation effects on 
particle production in pA collisions. At the same time, it will serve as a benchmark for 
assessing the hot medium effects on the quarkonium production in AA collisions. 

In the CGC framework, the quark-pair production cross section in pA collisions is 
obtained by Blaizot, Gelis and Venugopalan [23, 24], where pA is treated as a "dilute- 
dense" system and the cross section is evaluated at leading order in the coupling constant 
as = g'^/4:7i and the color charge density Pp in the proton, but in full orders with respect to 
the color charge density g'^pA = 0{1) in the nucleus. Adopting this formula, we previously 
evaluated the heavy quark production cross-sections in high-energy pA collisions to reveal 
their general features[10, 11]. 

Of particular importance is the x dependence of the unintegrated gluon distribution 
(uGD) in the nuclear target. We describe the x dependence of uGD with the Balitsky- 
Kovchegov (BK) equation including the running coupling corrections (rcBK equation) [25]. 
The BK equation can be derived from a more general equation of functional renormal- 
ization group by taking the mean-field approximation. The authors of [26, 27] analyzed 
the global data of deep inelastic scatterings (DIS) of e-|-p at HERA using the rcBK equa- 
tion, to obtain a set of constrained uGD of the proton at small x < 0.01. The resultant 
uGD has been applied to compute the particle production at the proton- (anti)proton 
colliders [28, 29]. We use this constrained uGD in this work. 

In forward particle production, the momentum fraction Xi of the gluons from the 
proton is not small, and the distribution may be better described with the ordinary 
collinear gluon distribution function. Accordingly, the expression of the cross section is 
obtained by taking the collinear limit fci^ — )■ for the gluon momentum from the proton 
in the hard matrix element. (This kind of asymmetric treatment is well known for the 
hadron production from CGC [30].) We will compare the quarkonium production in the 
collinear approximation on the proton side with the original one that involves the ki±_ 
dependent uGD for the proton. 

This paper is organized as follows: In Sec. 2 we review the quark pair production 



formula in the CGC framework and explain its collinear limit on the proton side. Then 
we incorporate the x-dependence of uGD using the rcBK evolution equation. In Sec. 3, we 
compute the quarkonium production cross section working in the color evaporation model, 
and show the transverse momentum spectra and nuclear modification factor for J/?/; at 
RHIC and LHC energies, and those for T at LHC energy. Rapidity and initial-saturation- 
scale dependences of the nuclear modification factor, and momentum broadening are also 
discussed. Sec. 4 is devoted to conclusion and outlook. 

2 Pair production cross-section 

In this section we summarize analytic expressions for the quark-pair production cross- 
section, derived in Ref. [24]. We also present its limit when the transverse momentum of 
the gluon coming from the proton is small, in order to make use of conventional collinear 
factorization for the proton. For more details, see Refs. [24, 11]. 

2.1 Pair cross-section in the large N limit 

In the CGC formalism, the proton-nucleus collision is described as a collision of two sets 
of color sources representing the large x degrees of freedom in the proton and the nucleus 
respectively. When they collide, these color sources produce a time-dependent classical 
color field, and this color field can in turn produce quark- ant iquark pairs. Both the 
classical color field and the quark-pair production amplitude can be calculated analytically 
when one of the projectiles is dilute and its density of color sources can be treated to the 
lowest order. This is the approximation we make in the description of proton-nucleus 
collisions in this framework, the proton naturally being the dilute projectile. 

For definiteness, let us denote by Pp and p^ the densities of color sources in the proton 
and the nucleus respectively, and that the proton moves in the +z direction and the 
nucleus in the —z direction. Then, the pair production amplitude reads [24] : 

(27r)2(27r)2 kf^ jax^ay^e e 

x^q){T,^{k,^,k^)[U{x^)eu\y^)] + Tg{k,^)[t'u'''{x^)]\v{p), (1) 
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^Thc momenta p and q of the produced particles have not been hsted among the arguments of these 
objects in order to make the equations more compact. 



with k2i_ = p^ + 9± ~ ^i±) and where C^{p + g, fci_|_) is the well-known Lipatov effective 
vertex^. The matrix U{x^) is the path-ordered exponential of the gauge field generated 
by the charge density p^ in the nucleus: 



U{x^) = V+ exp 



hOO 



i/ / dz^—YP^{z+,x±) -t 



(3) 



with t"" the SU(A^) generators in the fundamental representation, and U{xj_) is the same 
but in the adjoint representation. The two terms in the curly bracket of Eq. (1) may be 
interpreted as the processes where a gluon from the proton emits the pair either before 
or after the collision with the nucleus. The eikonal phases U and U describe the multiple 
scatterings of a quark or a gluon in the target nucleus. An important property of this 
amplitude is that the sum of the two terms in the bracket vanishes when ki± — )■ in 
agreement with Ward identities [24] - this property is essential in order to recover the 
limit of coUinear factorization on the proton side, which should hold since there is only a 
single scattering on the proton. 

The probability of pair production is then obtained by squaring the above amplitude, 
and by averaging the squared amplitude with the classical charge distributions of the pro- 
ton and the nucleus, lVp[x,pp] and VT^fx, p^]. The resulting expression generally involves 
the nuclear multi-parton correlators, e.g., {tr(^U{x±)t"'W{yj_)U{y'^)t"'W{x'j_)))y, where 
{■■■)y indicates that the average over p^ is performed with the distribution PV^[x,p4] 
evolved to the rapidity Y = ln(l/x). 2- and 3-point correlators are also needed in the 
calculation of the pair production probability, but they can be obtained as special cases 
of the 4-point function thanks to the identity U{x±)t"'U^{x±) = t'^U^"'{x±). This formula 
also leads to a general sum rule for the Fourier transforms of these correlators d'^^''^^, d'^.^'^ 
and 05/^ [24,11]. 

All these correlators can be evaluated in a closed form when the distribution of color 
sources, VF^[x, p^], is a Gaussian^ (see [24]), but the 4-point correlator has a very com- 
plicated expression which is quite hard to evaluate numerically. However, in the large A^ 
limit, it simplifies into 

~ _ _ _ /V2 

tT{Uix^)t^U\y^)U{u^)t^U\v^)) = —SAx^,v^)SAu±,y^) (4) 



with 



SA^^, y±) = ^tr(f/(x^)f/t(y^))^ . (5) 



^Its components are : 

C+{q,k,^) ^ ^^+q+ ; C;(g,fcix) ^ % -g- ; q(g,fci±) ^ -2kl+q^ . 
q ^ q+ ^ 

■^This distribution is a Gaussian in the McLerran-Vcnugopalan model for a large nucleus, and also 
in the asymptotically small x regime at very high energy. We therefore expect that this is a reasonable 



approximation for our work. 



As one can see, it is possible in this limit to write the 4-point (and also the 3-point) func- 
tion in terms of the 2-point function only, which simplifies considerably all the numerical 
calculations. Together with the translational invariance in the transverse plane, this fact 
makes the relation between 61^''^^ and ch^^'^ trivial. 

^ A,Y ^ A,Y 

By exploiting these relations between the 4-, 3- and 2-point correlators, the pair pro- 
duction probability at the impact parameter b, in the large N limit, can be written as [11] 
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(6) 



where we denote J^ = J d'^k±/(2'K)'^, d^ = N"^ — 1 the dimension of the adjoint repre- 
sentation of SU{N), and ki± = p± + qj_ — k2±. The variables yi^2 are the rapidities of 
the gluons that come from the proton and from the nucleus respectively. A shorthand 
notation for the squared matrix element is introduced as^ 



S(A;ij_,fc2±,fcj 



trd 

+trd 

+trd 



(^+m)T,,-(x/-m)7°Tj7° + h.c. 



(^+m)T,(2/-m)7°TV 
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In Eq. (6), '^p^yi is the uGD in the proton, and dc^^'^^^/'^^^J- is expressed in terms of the 
Fourier transform of the 3-point nuclear correlator as^ : 
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large A^ 4^^ 



Syiki) Sy{l^ -fcj 



"^In general, the Fourier transform of the 4-point function depends on three momentum variables: fc2_L, 
fcx and A;'j_ (see [24]). However, in the large N limit, this 4-point function is in fact given by (see [11]) 
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which allowed us to equate fe^ and h ^ in the squared amplitude and to perform directly the fej^ integra- 
tion. 

^In Eq. (6), this object appears in differential form with respect to the transverse coordinate X\_ 
because we are considering here the pair production probability at a fixed impact parameter h. When 
we integrate it over h in order to obtain the cross-section, we will have to integrate this correlator over 
the transverse area of the nucleus. 



where 5"^ (fc^) is the Fourier transform of Sy{x±). The X± dependence is rather weak for 
a large nucleus and may be treated implicitly through the variations of saturation scale 
Q^^(Xx) with X±. In the second line of Eq. (8), we have ignored the X^-dependence 
of 0^^'^ when doing the X± integration because the proton radius is small compared with 
that of a heavy nucleus {Rp <^ Ra)- ^^ thus have a compact expression for the quark 
production probability, but our formula involves a nuclear 3-point function 4>'^^y, which 
violates the usual form for /c^-factorization even in the leading order approximation [24, 
10]. 

The pair cross-section in the minimum-bias pA collision is obtained by integrating 
Eq. (6) over the impact parameter b. Dividing the cross-section with the total inelastic 
cross-section cr^^^^, which we estimate as cr^^^^ = 7r(i?^ + i?p)^ ~ Tri?^, we have the average 
multiplicity per event : 

dNqq _ 1 a^^N 1 /■ S(fclJ_^fc2±^fcjJ -g 

(9) 
Here we have introduced the 3-point function integrated over nuclear transverse area 

0ff(i±,fc±) = TiR\^S,{k^)S,{l^-k^), (10) 

which is related to the uGD (2-point function) of the nucleus by 0^'^(Z±) = /j^, 0^^^ (^±, k±). 
The proton uGD ^pp^y may be estimated by replacing the transverse area t^R\ and the 
amplitude Sy with those for the proton. 

2.2 Collinear limit on the proton side 

When the momentum fraction Xi probed in the proton is not small (e.g., Xi > 10"'^), 
and even more so in the forward rapidity region where Xi = 0{1), the typical transverse 
momentum of the gluons in the proton is much smaller than the transverse mass of the 
produced quark or antiquark, rrip^ ^ ki± = 0{A^^^). We can therefore neglect ki^ in 
the matrix element S in Eq. (7), and take the collinear approximation on the proton side. 
This limit is well defined thanks to the fact that the expression on the second line in the 
amplitude in Eq. (1) goes to zero as ki± -^ [24] : 

MM,P) = A-k^^ + Oikl^). (11) 

Thus, the amplitude squared S is quadratic in ki± when ki± — > 0, which cancels the factor 
kf_j_ in the denominator of Eq. (9). Note that the "vector" A in this formula contains 
spinors and Dirac matrices. In this approximation, we can write the integral in Eq. (9) as 

(12) 



where it is now implicit that ki± should not exceed the typical transverse momentum 
scale set by the produced final state. Using (Pki^ = ^d9id{kf^) and 
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we obtain 
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where we have now A;2± = Pj_ + Q±, and where we denote Scoii(fc2±, k±) = itrd(A^). The 
squared matrix element Hcou in the collinear approximation can be obtained by expanding 
the amplitude in Eq. (1) to linear order in the transverse momentum ki±: 



--coll --coll "T "coll "T "coll ' 



(15) 



with 



"coll 



"coll 



8p+g^ 



■"coll 



(p+ + g+)2(a^ + m?Y 

16 
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2 , M^+(?+)^2 
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-a. 



m -\ ^-r- — 7T^^a,± ■ {p^qj_ - q p±) 



(p + q) 
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In these formulas, we denote a± = q^^ — k±, and the squared invariant mass of the pair, 
{p + qY, is given by 



(p + qf = {p^ + q^) 



'p\ + m^ q\ + im? 
p+ q+ 



iP± + Qa 
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2.3 Correlators and energy evolution 

The dense nuclear distribution of gluons is encoded in the 3-point correlator (p'^f'-^, that 
appears in the pair cross-section (9) or (14). In the large N limit, it can be expressed in 
terms of the 2-point correlation 5'y(A;^), i.e. the forward scattering amplitude of a dipole 
in the fundamental representation. 

In the quasi-classical MV model, 5'y(fc) is independent of the rapidity variable Y, and 
it only includes the effects of the multiple scatterings of the quark-antiquark pair passing 
though the target nucleus, in the eikonal approximation. At large k±':^ Qg, the effect of 
multiple scatterings becomes small and leading twist results are recovered. 

The energy (rapidity) dependence of S'y (A;) arises via quantum fluctuations. Generally, 
the evolution equation for the 2-point function involves an inflnite hierarchy of multi-point 



correlators. However, it is known that in the hmit of a large number of colors A^ and of 
a large nucleus the energy evolution of Sy{k_i_) can be described by a closed mean-field 
equation known as the Balitsky-Kovchegov (BK) equation [31, 32]. This equation is an 
integro-differential equation that reads^ 



_d_ 
dY 



SAr±] 



dr^ilC{rjL,ri±) 



SArj 



-'S'y(ri_L)S'y(r2_L) 



where r± = ri± + r2± and }C{r±, ri±) is the evolution kernel (see below). Thus, with an 
appropriate initial condition at a certain a; = Xq, we can consistently treat the rapidity 
dependence of the cross-section by substituting into Eq. (8) the solution S'y(fc_|_) of the 
BK equation for x < Xq. 

It is also well known that the BK equation with a fixed coupling constant requires a 
very low value of as in order for the evolution of the saturation scale to be compatible with 
what one infers from HERA data, i.e., Ql(Y) ~ exp{XY) with A ~ 0.3 [6, 7]. It was argued 
that the next-to-leading order corrections to the BK equation would give the correct 
evolution speed with more reasonable values of a^ [33] . Indeed, it has been demonstrated 
recently in Refs. [25, 34] that the BK equation including the running coupling corrections 
in the kernel in Balitsky's prescription [35]: 



}Cir±,ri_L] 



27r2 



aJr 



a An 



1 + 



^2^2 



To 



a., r; 



ajr 



(19) 



makes the saturation scale behave compatible with HERA data, and the x-evolution 
equation becomes now a very useful tool (called rcBK equation) for phenomenology. 

Global fitting of the compiled HERA e+p data at x < Xq = 0.01 was performed with 
the rcBK equation in [26, 27]. Following their approach, we choose the initial condition 
at Yq = ln(l/xo) as 



5.0 (rx) 
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and the parameter values are listed in Table 1[29]. As for the running coupling constant 
in the evolution kernel we adopt the following form in the coordinate space: 



as{r') 



bnln 



4C^ 

^2^2 



1 -1 



(21) 



A constant a is introduced so as to freeze the coupling constant 
= afr. The non-Gaussian value 7 > 1 is preferred by the fitting. It is 



with 60 = 9/(47r) 

smoothly at as(oo 

argued in [36] that the value 7 > 1 suggests a possible importance of higher-order color 

correlations in the proton and is valid for moderate values of transverse momentum k_j_. 

We also list the McLerran-Venugopalan model 7 = 1 for comparison. 



^We have written it here in the approximation where the nucleus is translation invariant in the 
transverse plane. This is a reasonable approximation for a large nucleus, since the edges of the nucleus 
- where this invariancc is broken - give a comparatively small contribution to the cross-section. 
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QloJGeY' 
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afr 


C 


gins 

MV 


0.1597 
0.2 


1.118 
1 


1.0 

0.5 


2.47 
1 



Table 1: Parameter values for the two-point correlator Sy- A = 0.241 GeV is fixed. 



For xq < a; < 1, we extrapolate the function 
Ansatz [37]: 



/,g9.9 with the following phenomenological 



m,9 

A,Y 



{l±,k 



-L 



A,Yo 



{l_L,kj 



1 



X 



X J 



0.15 



(22) 



1 -a;o. 

In this formula, the power 4 for the factor 1 — x comes from the behavior at large x of 
the gluon distributions, as inferred from sum rules. Note that this extrapolation implies 
that the saturation scale is frozen at large x, which may lead to a harder fcx-spectrum 
for X > Xq than expected, possibly overestimating the Cronin peak. 

The saturation scale for a heavy nucleus at moderate values of x will be enhanced by 
a factor of the nuclear thickness TA^b) : Qlj^{x,b) oc TA^b) Qlp{x) [2]. As we consider 
only mean bias events in this work, we assume a simpler relation 



QIaM 



A'/-'Q 



s,p(^o) 



(23) 



where A is the mass number of the nucleus under consideration. This relation is assumed 
to be valid at the Xq where the initial condition for the rcBK equation is set. The rcBK 
equation then controls how the saturation scale of the nucleus evolves to lower values of x. 
Taking into account the uncertainty of the saturation scale for a heavy nucleus A ~ 200 
at X = Xq, we will consider the values in the range Q^^ = (4 — 6) x Ql^ at Xq = 0.01. 

In Fig. 1, we show the profile of uGDs, (Pp^y{k±) and (pA,y{k±), of the proton and the 
nucleus, respectively, at several values oi y = ln(xo/x), obtained by solving the rcBK 
equation with set glll8 in Table 1. We also show the result of the BK evolution with 
the fixed coupling constant a^ = 0.1, for comparison. This small coupling is necessary to 
keep the evolution speed compatible with the empirical value in the BK equation. We see 
that with increasing rapidity y the peak position (i.e., the saturation scale) drifts slightly 
faster in the rcBK evolution than that in the BK evolution with a^ = 0.1, and that the 
former yields a steeper /c^-slope than the latter. Note that the uGD should behave as 
l//c^ at large k± as is known in the LO double log approximation for BFKL or DGLAP 
equation, while the BFKL (equivalently BK in the linear regime) evolution gives harder 
k± spectrum. A dip structure around k^ = 3 GeV at x = a;o is caused by the parameter 
7 > 1, but is soon smeared out in the evolution. In the nucleus case, we assumed the 
initial saturation scale Q^q.a = ^Q'sQ,p at x = xq = 0.01. The uGD is more suppressed in 
low k± region and the peak position locates at larger k± than in the proton case, refiecting 
the stronger multiple scatterings in the nuclear target. 

We consider here the Xi^2 coverage of the charm pair production in the plane of the 
rapidity y and the transverse momentum P± of the pair at collision energies y^=200 
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Figure 1: Left: Evolution of the proton uGD, (l)y{k\)/{7TR'^N/4:as) obtained by solving 
the rcBK equation (solid) and the BK equation with fixed a^ = 0.1 (dashed) with the 
initial condition of set glll8. Right: the same with replacing Q^^ = GQ^Qp. 

GeV and 5.02 TeV in Fig. 2. Here we fix the pair's invariant mass M = 3.1 GeV, and 
draw the curves determined by Xi_2 = g^^{^/P± + M'^/y/s), on which either xi or X2 is 
constant. The kinematically disallowed region where xi^2 > 1 is indicated by the shaded 
area. We see that, at the RHIC energy, J/ip is produced from the gluons of moderate 
Xi^2 ~ 0.01 — 0.05 at mid-rapidities, while at forward rapidities ?/ ~ 2 the process gets 
sensitivity to the gluons at small X2 < 0.01. At the LHC energy, on the other hand, J/ip 
production is already sensitive to the small X2 gluon even at mid-rapidity, and at forward 
rapidity it probes X2 as low as ~ 10"'-^"^^ 

In the next section, we study the quarkonium production in GEM applied for the 
heavy-quark production cross-section (9) or (14) with the x-evolved uGD 0^*y(fc). GEM 
has been successful in describing the J/ip production in high-energy proton-proton (pp) 
collisions. 

3 Quarkonium production 

We estimate the J/tp production from the quark-pair production cross section within 

GEM[15]: 



da 



j/^ 



(PP^dy 



'J/^ 



AM% 



dM^ 



da. 



irrii 



d^P^dM^dy 



(24) 



where mc {Mn) is the charm quark (Z)-meson) mass. A phenomenological constant -Fj/^ 
represents the nonperturbative transition rate for the charm pairs, produced in the in- 
variant mass range from 2mc to the threshold 2Mo, to bound into a quarkonium. Its 
empirical value is around Fj/^ = 0.01-0.05[38]. Use of GEM for pA collisions assumes 
that the bound state formation occurs outside the target nucleus. We fix the threshold 
with Md (Mb) = 1.864 (5.280) GeV for 3/tp (T). 
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Figure 2: Kinematical coverage of the pair production in the plane of rapidity y and 
transverse momentum Pi for invariant mass M = 3.1 GeV at (a) y^=200 GeV and (b) 
^/s=5.02 TeV. Shown are the curves of constant Xi.2 = (a/-PJ + M^/-y/i)e^^. The shaded 
region is kinematically forbidden. 

A remark is here in order. In the pair-production cross section (9) the inelastic cross 
section estimated as 7rR\ in the denominator effectively cancels out with the same factor 
in (p'^^y, and that the cross section is now proportional to the effective transverse area 
TrRp of the proton appearing in ^Pp^y In the following calculations, we choose the proton 
size Rp = 0.9 fm and the J/ip formation fraction -Fj/^ = 0.02 as representative values. 
One should keep in mind that the absolute normalization of the cross section depends on 
these parameters. We also cancels a^ in front of the cross section by a^ appearing in the 
denominator in (f)A,y and ^Pp^y In the case of coUinear approximation on the proton side, 
we set as = 0.2 in this paper. 

3.1 Transverse momentum spectrum of J/i/j 

RHIC 

We first show in Fig. 3 the transverse momentum spectrum of the produced J/ip in pp 
collisions ^ at ^/s = 200 GeV, using the uGD set glll8 given in Table 1. The upper (lower) 
curve of each band indicates the result with charm quark mass rric = 1.2 (1.5) GeV. In the 
collinear approximation on the larger-xi side, we adopt CTEQ6L0 parametrization[40], 
and the band in Fig. 3 includes the change of the factorization scale from 2M± to M±/2 
with M± = (M^ + P^)^/^, where M is the pair's invariant mass. 

As mentioned above, the quarkonium production at mid-rapidity \y\ < 0.35 is largely 
determined by the gluon distributions at moderate Xi 2 ^ 0.01. Then, we notice a difficulty 
with set glll8: the peculiar dip structure of glll8 seen in Fig. 1 remains in the J/ip 
spectrum as a similar dip around P^ ~ 2 GeV, which must be an artifact of this initial 



^Strictly speaking, treating pp (at mid-rapidity) as a dilute-dense system is not legitimate but we need 
the pp cross sections for studying the so-called nuclear modification factor RpA (see text). 
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Figure 3: Transverse momentum spectrum of J/ip in di-lepton channel in pp collisions 
at ^ = 200 GeV for rapidity ranges (a) \y\ < 0.35 and (b) 1.2 < y < 2.2. GEM 
model results using the pair production (9) with sets MV and glll8 are shown in gray 
and doubly-hatched bands, respectively, and the result using coUinear approximation (14) 
with set glll8 is in hatched band. The upper (lower) curve of the band corresponds to 
the result with rric = 1.2 (1.5) GeV, and the scale of pdf is chosen at 2M± {M±/2) in the 
coUinear approximation. Data from [39]. 



condition. In contrast, we don't see such a structure with the MV initial condition. At 
forward-rapidity 1.2 < y < 2.2, the dip is smeared to be less noticeable by the imbalance 
between xi and X2 and by the X2 evolution of the uGD. As a whole, the P± spectrum 
obtained with set glll8 is closer to the observed data [39] than with set MV. In this pp 
case, the coUinear approximation on the large-Xi side does not improve the description of 
the data. The k± kick from only the one of the protons cannot give enough Pi for the 
pair. 

In Fig. 4 shown is the transverse momentum spectrum of the J/ip in pA collisions 
in our model. We set the initial saturation scale of the uGD for the heavy nucleus as 
Q'so A — 6Qso,p at X = ajQ. The upper (lower) curve of each band indicates the result with 
me = 1.2 (1.5) GeV. We overlay d-Au data observed by PHENIX at ^s = 200 GeV[41], 
presuming here that the difference between pA and dA results only in normalization 
difference of order 0(1) ^. We find that Pi-dependence of J/ip production is better 
described with set glllS ^ than that with set MV. Indeed, here the coUinear approximation 
on the proton side apparently gives a better description of the data both at mid- and 
forward-rapidity regions. However, at forward rapidities, where we approach the small-X2 
region and the kinematical boundary for xi at the same time (see Fig. 2), we expect a 
nontrivial interplay between large xi and small X2- Besides the saturation dynamics of X2 
gluons, one may need to consider other physics such as energy loss of large-xi gluons in 



^Recall that our model already has an uncertainty of 0(1) in the normalization of the uGD. 
^Possible dip structure from the proton uGD is smeared out here in Fig. 4 by the multiple scattering 
effects in the nuclear uGD. 
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(a) Vs = 200 GeV, |y| < 0.35 



(b) Vs = 200 GeV, 1 .2 < y < 2.2 
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Figure 4: Transverse momentuin spectrum of ^/ip in di-lepton channel in pA collisions at 
\/s = 200 GeV for rapidity ranges (a) \y\ < 0.35 and (b) 1.2 < y < 2.2. Notations are the 
same as in Fig. 3. Data in d-Au collisions[41] are overlaid for comparison. 



the heavy target [43] in order to understand the P± spectrum of J ftp in the very forward 
region. These effects are not included in our present treatment. 

We notice in Fig. 4 (b) that the J/ip production is more suppressed nearly by one 
order of magnitude in the collinear approximation than those in the full calculation. This 
is caused by a difference in the large xi behavior of gluon distributions on the proton side. 
As X — ;• 1, the CTEQ gluon distribution decreases much more rapidly than our model 
uGD v^p,j/, which is assumed as oc (1 — x)^. Furthermore, in the collinear approximation, 
the pair's P± is entirely provided from the nucleus side, P± = k2, and uGD (pA^y for the 
heavy target is more suppressed at low /c2 by multiple scatterings. 

Now let us take a ratio of the cross section of J/ip in pA collisions to that in pp colli- 
sions, which is called nuclear modification factor R^a- We expect that model uncertainties 
cancel out to some extent in the ratio. We define -RpA for J/ip in our model as 



R 



da.y^,/(PP±dy\ 



pA 



pA 



A^cou daj/^/d'^P±dy\ 



pp 



(25) 



A^'^ because the 



Here we set the number of nucleon-nucleon collisions in pA to A'^coii 
uGD (l)A,yoikj_) scales as (Qso)^ ^ ^'^^^ ^^ large k±. 

In Fig. 5 we compare the model results for RpA at ^/s = 200 GeV with the data of 
-RdAu- Note that the projectile is different between the model calculation and the data. 
The notations are the same as in Fig. 3. We stress here that RpA is indeed little dependent 
on the choice of the quark mass and factorization scale. Unfortunately, however, one 
immediately recognizes an unphysically strong Cronin peak in the model calculations 
with set glll8 both at mid- and forward rapidities, which is obviously caused by the 
dip seen in the pp collisions (Fig. 3). In contrast, the RpA result with set MV looks 
more reasonable; we see a moderate Cronin peak at mid-rapidity due to the multiple 
scatterings, while it almost disappears at forward rapidity y ~ 2 by the X2 evolution. In 
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(a) Vs = 200 GeV, |y| < 0.35 



(b) Vs = 200 GeV, 1 .2 < y < 2.2 





Figure 5: The ratio of J/^ productions in pA and pp collisions -RpA(-Pi) at ^/s = 200 
GeV for (a) \y\ < 0.35 and (b) 1.2 < y < 2.2. The results with uGD sets MV and glll8 
are shown in gray and doubly-hatched bands, respectively, and the result in collinear 
approximation with set glll8 is shown in a hatched band. Notations are the same as in 
Fig. 3. Data of RdAu taken from [41]. 



low- Pi region, we also notice a too strong suppression than the experimental data. This 
would imply the importance of the fragmentation process in the formation of J ftp, which 
is missing in a simple GEM treatment. 

To summarize the results at RHIG energy, the J/ip production spectrum is sensitive 
to the moderate value of Xi^2; where the initial condition for the x-evolution is set. We 
have a difficulty to describe the pp data and therefore the ratio RpA with the uGD glll8 
constrained at x < 0.01. In contrast the set MV gives more reasonable behavior for 
-RpA. In pA collisions the P± spectrum is better described with set glllS at mid- and 
forward-rapidities. In forward rapidity, the observed P± slope is still steeper than the 
model, hinting other effects such as a possible energy loss of the large-a;i gluon from the 
proton. Actually i?dA of J/t/.' at RHIG energy has been studied in several approaches 
(e.g.) with introducing nuclear parton distribution and nuclear absorption effects to a 
J/ijj production model for pp[42], or with taking account of the multiple scatterings and 
energy loss of the projectile gluons[43]. 

LHC 

Now we compute the J ftp production at the LHG energy, where we expect that a wider 
X2-evolution of uGD on the nucleus side will manifest in the quarkonium spectrum. In 
fact, both xi^2 are small (~ 10"^ < xq) already in mid-rapidity production of the charm 
pair as seen in Fig. 2, and as moving to larger rapidities we can probe smaller values of 
X2 on the nucleus side down to X2 ~ 10~^. 

We show in Fig. 6 the J/ip cross section in pp collision at -y/s = 7 TeV, obtained 
in GEM from charm quark spectrum (9). Notations are the same as in the case of the 
RHIG energy. In order to assess the uncertainty, we again vary the charm quark mass from 
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(a) Vs = 7 TeV, ly| < 0.9 



(b) Vs = 7 TeV, 2.5 < y < 4 
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Figure 6: Differential J/'?/' yield in pp collisions at y/s = 7 TeV for (a) \y\ < 0.9 and (b) 
2.5 < ?/ < 4. Notations are the same as in Fig. 3. Data from [44]. 



(a) Vs = 5.02 TeV, -1 .365 < y < 0.435 



(b) Vs = 5.02 TeV, 2.035 < y < 3.535 
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Figure 7: Transverse momentum spectrum of ^/ip in di-lepton channel in pA collisions at 
i/s = 5.02 TeV for (a) —1.4 < \y\ < 0.4 and (b) 2 < y < 3.5. Notations are the same as 
in Fig. 3. 
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(a) <s = 5.02 TeV, -1 .365 < y < 0.435 



(b) Vs = 5.02 TeV, 2.035 < y < 3.535 
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Figure 8: The ratio RpA{P±) for J/tp at y/s = 5.02 TeV for (a) 
2 < y < 3.5. Notations are the same as in Fig. 3. 



-1.4 < y < OA and (b) 



rric = 1.2 to 1.5 GeV, and change the factorization scale from 2M_|_ to M^/2 in the colhnear 
approximation. The observed data[44] is fairly well reproduced with set glll8 in this P^ 
region both at \y\ < 0.9 and 2.5 < y < 4, indicating that y-dependence is appropriately 
captured by x evolution of uGD. The Pi slope in the colhnear approximation (14) with 
set glllS seems to be slightly off the data, while the full result with set MV gives harder 
Pi spectrum. The situation is expected to be similar in pp collisions at ^/s = 5.02 TeV. 

Results in pA collisions at ^/s = 5.02 TeV are plotted at mid- and forward-rapidities 
in Fig. 7. The MV initial condition gives a harder spectrum of J/ip than glll8. But their 
Pi slopes become almost the same at Pi > 10 GeV, hinting the same BFKL tail of uGD 
generated during the evolution. Compared to the case at the lower energy ^/s = 200 GeV, 
the colhnear approximation (with set glllS) results in the spectral shape rather similar 
to the full result at this energy ^/s = 5.02 TeV, where the colhnear approximation on 
the proton side would be more suitable since the saturation scale of the nucleus is much 
larger than that of the proton: Qg^(x2) ^ Qlp{xi), especially in the forward region. 

We show in Fig. 8 the ratio PpA of J/ip as a function of Pi at y/s = 5.02 TeV. We 
have assumed iVcou = A'^^^ as mentioned before. We find that each band almost collapses 
into a single line, which means that the ratio PpA is insensitive to the variation of the 
charm quark mass (and the factorization scale in the colhnear approximation) within the 
range considered here. 

At mid-rapidities (Fig. 8 (a)), we see that the ratio PpA of J/ip production is suppressed 
at low Pi, while it approaches unity at higher Pi for both sets of glll8 and MV. In 
the colhnear approximation on the proton side, PpA shows a Cronin-like peak around 
Pi ~ 4 GeV and remains larger than unity at larger Pi, which largely reflects "PpA for 
4>A,y" at the gluon level. At forward rapidities (Fig. 8 (b)), however, this difference due 
to different uGD sets and approximations becomes much weaker to yield a systematic 
suppression as a function of Pi for all three cases. 

We examine the initial-scale (Q^o a) dependence of the ratio PpA in Fig. 9, by plotting 
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a) <s = 5.02 TeV, -1 .365 < y < 0.435 




(b) \'s = 5.02 TeV, 2.035 < y < 3.535 
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Figure 9: Initial-scale dependence of the ratio RpA{P±) for J/tp at (a) mid- and (b) 
forward-rapidities at i/i = 5.02 TeV. Q^q.a i^ set to 4(5so,p (upper) and 6Q%p (lower). 

the results with the saturation scale Q^q.a — ^Q'so,p (upper) and GQ^Qp (lower) in Eq. (20) 
at a; = Xq. It is found that the Q^q ^ dependence of RpA is relatively weak within this range. 
At low Px we have strong suppression, but one should keep in mind that this suppression 
may be filled to some extent by the nonperturbative fragmentation in forming J/ip, as is 
inferred from the discussion on Fig. 5. 

To summarize the result at LHC energy, we can probe here a wide X2-evolution of the 
uGD 0^,1/2(^2) through the J/ip production, and the ratio -RpA will be a good indicator 
for it. 

3.2 T production at the LHC 

Next we consider T(IS') production. Non-linear effects are generally suppressed by the 
inverse power of the heavy quark mass. However, since the bottom quark mass m?, is just 
three times as heavy as the charm quark mass rric, the relevant value of x for the T(15') 
production becomes larger by the same factor at low Pj_, as compared to the J ftp. At the 
LHC energy, this x value may be still small enough for multiple scatterings and saturation 
to be relevant in the T production. 

We plot the P± spectrum of T(IS') in pp and pA collisions at ^/s = 7 and 5.02 
TeV, respectively, in Figs. 10 and 11, together with the data measured by ATLAS and 
LHCb[45, 46] for the pp case. Here we have chosen the CEM parameter as -Ft(is) = 0.01, 
and varied mi, from 4.5 to 4.8 GeV. Other notations are the same as in the J/ip case. In 
pp collisions, the coincidence between the model and the data for T(IS') state is not as 
good as that for J ftp at low P± and at forward rapidity. 

We present in Fig. 12 the nuclear modification factor Rp^ for T(15') as a function of 
P±. The model uncertainty from the quark mass value and the factorization scale would 
cancel out by taking the ratio of the cross-sections in the pp and pA collisions. Indeed, 
each band collapses into a thin line whose width is almost unnoticeable. 
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Vs = 7 TeV, ly| < 1 .2 
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Figure 10: Transverse momentum spectrum of T(15') in di-lepton channel in pp collisions 
at ^/s = 7 TeV for (a) \y\ < 1.2 and (b) 4 < y < 4.5. Notations are the same as in Fig. 3. 
Data from [45, 46]. 
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(b) ^'s = 5.02 TsV, 2.035 < y < 3.535 
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Figure 11: Transverse momentum spectrum of T(15') in di-lepton channel in pA collisions 
at ^/s = 5.02 TeV for (a) —lA<y< 0.4 and (h) 2 < y < 3.5. Notations are the same as 
in Fig. 3. 
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(a) <s = 5.02 TeV, -1 .365 < y < 0.435 



(b) Vs = 5.02 TeV, 2.035 < y < 3.535 
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Figure 12: The ratio i?pA for T(IS') at ^/s = 5.02 TeV as a function of the transverse 
momentum, (a) — 1.4 < y < 0.4 and (h) 2 < y < 3.5. Notations are the same as in Fig. 3. 

This result for T(IS') is quahtatively very similar to that for J/tp. At mid-rapidity, 
we see a suppression Rp\ in low Pj_ region below 5 GeV, while it turns back to unity 
at larger P^. Only in the coUinear approximation, we see the Cronin-like enhancement, 
which is largely caused by the dip structure in the proton uGD at moderate xi. At 
forward rapidities 2 < y < 3.5, the T production is suppressed in a wide P± region from 
to 20 GeV, irrespective of the model uGD's, glll8 or MV, or of the use of collinear 
approximation. In the forward region, T(IS') production has the sensitivity to the small-x 
evolution of uGD in the nucleus. 

We have also checked the initial-scale (Q^oa) dependence of RpA for T(IS') by com- 
paring the result with Q^q.a = ^Q'so,p ^^d 6(5so,p to find that the change is very similar to 
the case with J/ip (Fig. 9). 

3.3 Rapidity dependence of RpA of J/i/j and T 

We study the rapidity dependence of the ratio RpA integrated over P^. The computation 
is performed with set glll8. In Fig. 13 shown is RpA{y) of J/ip at -y/s = 0.2 and 5.02 
TeV, together with that of T(IS') for the latter. Note that our assumption of dilute-dense 
colhding system applies only in the positive rapidity region {y > 0), especially for pp, 
which is needed in the denominator of PpA- 

We see systematically a stronger suppression of PpA as the rapidity increases both 
at RHIC and LHC energies. This is in accord with x-evolution of uGD in the heavy 
target. R^Aiy) of J/V' flattens out at y < 1 at RHIC energy because the J/-^ is produced 
there by the gluons with X2 > Xq and we freeze the saturation scale to its initial value 
Qlix > xO) = Ql,. 

Comparing the results of J/ip and T(IS') at LHC, we note that the suppression of 
T(15') is smaller than that of J/ip, but is still significant to be observed. It would be 
quite important to study these systematics in experimental data in order to quantify the 
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Figure 13: Nuclear modification factor Rpx for J/^ integrated over P± as a function of 
rapidity, in pA collisions at (a) a/s = 200 GeV and (b) y/s = 5.02 TeV. The band includes 
uncertainty for m^ = 1.2 GeV to 1.5 GeV and Q^^^^ = (4 _ 6)Ql^^^. In (b), RpA of T{1S) 
is also shown where we vary rrih = 4.5 to 4.8 GeV. RHIC data from [47]. 



saturation effects in the heavy nuclear target. 



3.4 QIqa dependence of Rp^ 

It would be interesting to study the dependence of -RpA on the saturation scale parameter 
(5so A) which may be translated to the effective thickness of the target. We compute RpA 
of J/ip and T(IS') integrated over P± as a function of Q^q ^ at several values of y. We fix 
here the uGD set glll8 and the quark masses as m^ = 1-5 GeV and mi, = 4.8 GeV. In 
Fig. 14 we plot RpA of J/ip at y/s = 0.2 and 5.02 TeV. We found that for each rapidity 
Q^Q ^-dependence of RpA can be fitted nicely by a model function: 



R 



pA 



(6 + Q2 



(26) 



sO,A) 



with a, b and a being parameters. This functional form is motivated by QCD analog of 
superpenetration of a electron-positron pair through a medium[48, 11]. We show T(IS') 
result in Fig. 15 with fitted curves. The stronger suppression at the larger value of Q^q ^ 
is naturally understood as a result of stronger multiple scatterings and saturation effects 
in the heavier target. 

Energy and rapidity dependences may be qualitatively inferred through the increase of 
Q'i A^y) with increasing y. Thus we tried to fit the rapidity dependence of -RpA by replacing 
in Eq. (26) Q^g A ~^ Q'sOA^^^ with a free parameter A, but it was only unsatisfactory. We 
remark here that quarkonium suppression due to parton saturation in our treatment is 
twofold: a relative depletion of the gluon source and multiple scatterings of the quark 
pair in the target. The latter disturbs the boundstate formation, by increasing the pair's 
invariant mass on average in CEM[49]. It appears hard to describe energy and rapidity 
dependence of the suppression at the same time through a single function Q 
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Figure 14: Nuclear modification factor RpA for J/ip as a function of Q% a at ?/ = 0, 1,2 
and 3 at ^/s = 200 GeV (left) and ^/s = 5.02 TeV (right). Fitted curves are also shown. 
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Figure 15: Nuclear modification factor RpA for T(IS') as a function of Q^ at y = 0, 1,2 
and 4 at Vs = 5.02 TeV. 
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Figure 16: Mean transverse momentum square A(P^)pA of J/^ as a function of Q^g A ^^ 

'l] is 



^s = 200 GeV (left) and ^s = 5.02 TeV (right). Fit with a form a[{Q%JQ 
also shown. 
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3.5 Pl broadening 

Finally, we study the mean transverse momentum of quarkonium in pA collisions. The 
momentum broadening in the nuclear target has been discussed in the literature [15]. In 
our framework, the multiple scatterings of the incident gluon and the produced quark 
pair in the nuclear target, encoded in U and U terms in Eq. (1) respectively, cause the 
momentum broadening of the pair. Typical momentum transfer of the multiple scatterings 
in the nucleus should be characterized by the saturation scale Qs,a{x2)- We define here 
the broadening of P± as the deviation of the mean transverse momentum squared (Pf) 
of J ftp in pA collisions from that in pp collisions: 



A(Pl)pA = (Pl)pA - (Pi) 



PP 



/ dapA f rfo-pp 



(27) 



In Fig. 16 we plot A{Pj_)pA as a function of Q^q.a- We use uGD set glll8 with the 
quark masses rric = 1.5 GeV and rub = 4.8 GeV. We have found that for each rapidity the 
Qlo A dependence of the broadening can be fitted in a simple form: 



A(Pl) 



pA 



^[iQlo,A/Qlo,p) 



1] 



(2^ 



Iq,, which 



with a and a being parameters. 

At y/s = 200 GeV, the broadening at mid-rapidity is obviously linear in Q 
indicates the random walk nature of the multiple scatterings in the momentum space. 
In the forward region, we naively expected an increase of the mean momentum by the 
stronger scatterings, but actually found the opposite, i.e., a decrease from the mid-rapidity 
value. We interpret this as the effect of kinematical boundary of xi in the forward region 
(see Fig. 2). 

The measured value of A(Pl)dA at RHIC[41] seems to be smaller by a factor of 5 than 
that in Fig. 16, if we naively translate Q^q ^ to the centrality parameter A'^cou evaluated for 
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Vs = 5.02 TeV, rrih = 4.8 GeV 




Figure 17: Mean transverse momentum square A(P^)pA of T(IS') as a function of Ql at 
v/i = 5.02TeV. Fit with a form a[((5^o,A/<5so,p)" - 1] is also shown. 

dAu colhsions. This strong broadening originates probably from the fact that our model 
has too hard P± spectrum at RHIC energy. But it is at least consistent with data that 
P± broadening at forward rapidities ~ 2 is weaker than that at mid-rapidity y ~ 0. 

At ^/s = 5.02 TeV, a wider phase space opens up and we instead see an increase of the 
mean momentum of J/ip as moving to the forward rapidity region. We have checked that 
A(P^)pA gets back to be smaller aX y = 6 than that of mid-rapidity, just as seen in the 
case of a/s = 200 GeV. Nonlinear dependence on Q^q ^ may imply the different evolution 
speed of multiple-scattering strength for different initial values Q^o a- The result for T(IS') 
in Fig. 17 is similar to the J/ip case, but interestingly the broadening becomes more 
remarkable; The heavier bottom quark pair can acquire the larger transverse momentum 
P± in multiple scatterings before going beyond the threshold set on the pair's invariant 
mass M^ < 4M| 

4 Conclusion and Outlook 

Quarkonium production in proton-lead collisions provides us with a good opportunity 
to study the saturation phenomenon in the incident nucleus thanks to the wide kinetic 
reach at the LHC We have computed the J/tp and T(15') production in pA collision 
at collider energies within GEM based on the GGG quark pair production, and have 
discussed sensitivity of the quarkonium observables to the parton saturation in the target 
nucleus. We have presented the calculations with the uGD set glll8 which is constrained 
with DIS data at x < Xq = 0.01, with and without the coUinear approximation, and the 
calculations with the model uGD set MV for comparison. 

At the RHIG energy ^/s = 200 GeV, the J/ip at mid-rapidity is produced not from 
the small-x gluons, but from the moderate-x gluons, and the Pi spectrum in pp collisions 
is unfortunately sensitive to a unphysical dip structure of the uGD set glll8, which was 
constrained only for x < xq. We need better extrapolation of our framework to x > xq. 
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In pA collisions, multiple scatterings smear out the dip of the uGD and the P±_ spectrum 
of ^/ip becomes closer to the observed one in dAu collisions. 

At the LHC energy y/s = 5.02 TeV, the small-x gluons dominate the charm production, 
and we have found that our model with the uGD set glll8 works for J/'0 production in pp 
collisions both at mid- and forward-rapidities. Then we have shown our model prediction 
on J ftp production in pA collisions. The ratio RpA{P±) for J ftp shows a suppression for 
-Pi ^ 5 GeV at mid-rapidity due to saturation effects, and it is further suppressed in 
wider range of Pi as moving to forward rapidities. 

We have also shown that the T(IS') production in pA collisions at the LHC has a good 
sensitivity to the gluon saturation of the nucleus, provided that the effect is smaller than 
that in the J/ip case. In our model, when integrated over P±, the ratio RpA{y) for T at 
the LHC shows a suppression similar to that of J/ip at RHIC energy. 

Transverse momentum broadening A(Pj) of the quarkonium shows an increasing be- 
havior as a function of QIq ^. Because our model gives harder P± spectrum than the data, 
the broadening is likely to be overestimated at RHIC energy in our calculation. However, 
it is still interesting to notice that at RHIC energy the broadening becomes weaker at 
forward rapidity than at mid-rapidity due to the kinematical constraint. At LHC en- 
ergy, on the other hand, we expect the increase of the broadening at forward rapidities 
because of the larger saturation scale Ql{x2) and wider kinematic coverage of the LHC. 
Transverse momentum broadening is also investigated recently by taking account of the 
multiple scatterings in the target in [50]. 

In conclusion, we have numerically studied the quarkonium production in pA collisions 
at the RHIC and LHC, within CEM based on the CGC quark-pair production energies, 
and have quantifies the effects of saturation and multiple scatterings in the target nucleus 
on the J/ip and T observables. Comparison of our results with experimental data at 
the LHC must be very important to access the relevance of the saturation physics in the 
quarkonium production. 

In this work we have employed CEM to describe the nonperturbative formation of the 
quarkonia. In fact, quarkonium formation is one of the challenges in QCD, even in pp 
collisions, and CEM replaces this just by a probability constant -Fj/^,. Non-relativistic 
QCD framework has recently been extended to NLO, which improves our understanding 
of quarkonium production significantly [51, 52]. As a first step in this direction we plan to 
match the quark-pair production from CGC onto the non-relativistic QCD approach. 

We can investigate also the production of open heavy flavor mesons in our framework. 
Modification of P± spectrum and DD correlations will also provide very useful information 
on the saturation in the target nucleus and a good benchmark for the energy loss and 
collective flow measurements of D and B mesons. We will report this elsewhere [53]. 
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